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Abstract. In this work we derive and analyze coarse-grained descriptions of self- 
propelled particles with selective attraction-repulsion interaction, where individuals 
may respond differently to their neighbours depending on their relative state of motion 
(approach versus movement away). Based on the formulation of a nonlinear Fokkcr- 
Planck equation, we derive a kinetic description of the system dynamics in terms of 
equations for the Fourier modes of a one-particle density function. This approach allows 
effective numerical investigation of the stability of possible solutions of the system. The 
detailed analysis of the interaction integrals entering the equations demonstrates that 
divergences at small wavelengths can appear at arbitrary expansion orders. 

Further on, we also derive a hydrodynamic theory by performing a closure at the 
level of the second Fourier mode of the one-particle density function. We show that 
the general form of equations is in agreement with the theory formulated by Toner and 
Tu. 

Finally, we compare our analytical predictions on the stability of the disordered 
homogeneous solution with results of individual-based simulations. They show good 
agreement for sufficiently large densities and non-negligible short-ranged repulsion. 
Disagreements of numerical results and the hydrodynamic theory for weak short-ranged 
repulsion reveal the existence of a previously unknown phase of the model consisting 
of dense, nematically aligned filaments, which cannot be accounted for by the present 
Toner and Tu type theory of polar active matter. 
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1. Introduction 

In recent decades, there has been an increased research focus on far-from-equilibrium 
systems in biology and physics which is referred to as "active matter". The relevant 
length scales of such systems span several orders of magnitude. They range from the 
(sub-)micrometer scale governing the dynamics of individual active units in motility 
assays in- vitro [I] and the actin cortex in- vivo [2], via the mesoscopic length scales 
of interest in collective dynamics of large bacterial ensembles [3]H] and artificial self- 
propelled particles [5], up to the macroscopic scales of driven granular matter [6l[7] or 
flocks of birds [8], schools of fish [9], and swarms of insects P^0|fTTj . where the spatial 
dimensions can be in the order of kilometers. 

Despite the apparent variety, all these systems share the fundamental property 
of local uptake and/or conversion of internal energy into kinetic energy of motion by 
its individual units. This - together with additional interactions between those units 
- distinguishes this systems from related equilibrium systems and yields fascinating 
examples of self-organization and collective dynamics. 

The question of universal properties of such active matter systems from the 
statistical physics point of view is a vibrant research field. To a large extent, it 
was initiated by the numerical study of a minimal, individual-based model of active 
matter published by Vicsek et al in 1995 [12]. Shortly after this publication, Toner 
and Tu made a seminal contribution by formulating the hydrodynamic equations of 
polar active matter at largest relevant length and time scales purely based on symmetry 
arguments [T31IH] . The analysis of this generic equations, as well as their counterparts 
for nematic order, improved our understanding of the fundamental properties of active 
matter, such as the existence of long range order or giant number fluctuations [T51IT6] . 

However, the direct derivation of a hydrodynamic theory of the Toner and Tu 
type from microscopic models of active matter was a long standing problem. Only 
recently, such a link between microscopic parameters determining the dynamics of 
individual active units and parameters governing the macroscopic flow of active 
matter was established by formulating kinetic equations for minimal models of self- 
propelled particles with velocity-alignment [TTl - [2~0] and self-propelled aligning rods [21] . 
Furthermore, coarse-grained descriptions for active particles with variable speeds and 
velocity-alignment were derived in [22H26] . 

In this paper, we will derive a kinetic and a hydrodynamic description for self- 
propelled agents interacting via a selective attraction and repulsion interaction. A 
corresponding model was recently introduced to describe the onset of collective motion 
in insect swarms and is directly motivated by response of individual agents to looming 
visual stimuli [271[28]. The model shows different phases including large scale collective 
motion, a disordered clustering phase and a nematic phase despite the absence of an 
explicit velocity-alignment interaction. 

First, we will introduce the microscopic, individual-based model in terms of 
stochastic differential equations. Then we will proceed with the discussion of a kinetic 
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description of the collective dynamics based on the nonlinear Fokker-Planck equation for 
the one-particle density function, which allows efficient numerical analysis of the stability 
of solutions of the Fokker-Planck equation (especially the spatially homogeneous, 
disordered state) in Fourier space. Further on, we will derive a hydrodynamic theory 
for self-propelled particles with selective attraction-repulsion interaction, which yields 
hydrodynamic equations in agreement with the theory by Toner and Tu. A direct 
comparison of the kinetic approach, which in principle can be considered up to arbitrary 
accuracy, to the hydrodynamic theory, which corresponds to a closure at the level of the 
second Fourier mode of the probability density function, reveals the range of validity 
of the hydrodynamic equations at large wavenumbers. Moreover, the analysis provides 
insights into the origin of unphysical divergences at large wavenumbers related to the 
approximations used namely the usage of Taylor polynomials. Finally, we compare the 
results of the kinetic and hydrodynamic theory with direct numerical simulations of the 
individual-based model. 



2. Microscopic Model 

We consider TV self-propelled particles of mass m = 1 in two spatial dimensions, so- 
called agents. Each individual moves at a constant speed s , thus the velocity vector of 
each agent is determined by its polar orientation angle ipi. The equations of motion for 
the positions r» and the polar orientation angles (pi read: 

dr< / cos ipi(t) 

— = s e iih (t) = s \ . , 2.1a 

at \ sm (pi(t) 

d(fi 1 



- [F^+ v^&W) • (2.16) 



dt s 

Here, F i)tp = F$ e ijV , is the projection of an effective social force vector F,i on the angular 
degree of freedom (ft, which induces a turning behaviour of the focal individual due to 
interaction with others. e i>ip = (— sin tp iy cos (pi) T is the angular unit vector perpendicular 

to Bi jh . 

The second term within the brackets in ( 12.1ft!) stands for random angular noise with 



intensity D v . £j(t) are independent, Gaussian random processes with vanishing mean 
and temporal ^-correlations, i.e. (£i(t)£j(t + r)) = 5ij S(t) (Gaussian white noise). 
The total social force is given by a sum of three components: 

fi,r ~\~ fi,m fi,a- (2-2) 

The first term represents a short-ranged repulsion allowing for finite sized agents. It 
reads 

N 

Ar = Hr(rji)rji6{l c - rji) (2.3) 

3=1 

with /J- r ( r ji) ^ being a repulsive turning rate, which in general depends on distance 
Tji = \rji\ = \vj — Vi\ between two particles. We assume that this repulsive interaction 
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is strictly short-ranged, i.e. it vanishes above a finite repulsive radius l c , as indicated by 
the Heaviside (unit-step) function 6(x). 
The other two forces read: 



fi,m = S ^^m{rji) \vji\ rji6{l s - rji)9(rji - l c )9(vji), (2.4a) 

5=1 
N 

fi,a = ^2 Vaira) \vji\ fji6(l s - rji)0(rji - l c )6(-Vji). (2.46) 

5=1 

These forces can be considered as a sum over pairwise interactions, which act always 
along the unit vector rji = (rj — Ti)jr,ji pointing towards the center of mass of the 
neighbouring particle j. The corresponding response strengths fi a ,m( r ji) are distance 
dependent and may in general be both positive (attraction) and negative (repulsion). 
Furthermore, the overall response to other individuals is assumed to vanish above a 
finite sensory range l s : fi a , m ( r ji > h) = 0, whereby l s > l c . 

The decisive factor in the distinction of the two forces is the sign of the relative 
velocity Vji defined by the temporal derivative of the distance between particles i 
and j, and equals the projection of the velocity difference Vjj = Vj — of neighbour 
j and the focal individual i on the relative position unit vector r^. fi^ a is the 
response to approaching individuals characterized by a negative relative velocity Vji < 0, 
whereas fi >rn is the corresponding response to moving away (or receding) individuals 
characterized by positive relative velocity Vji > 0. Both force terms are proportional 
to the absolute value of the relative velocity, leading to stronger responses to faster 
approaching or receding individuals. Using the definition Vji = fji, one obtains the 
relative velocity 



Vji{ip u cpj, atji) = 2s sin ^ ^ — - acjA sin ( — - ¥± j (2.5) 

as a function of the velocity angles <fii, <pj and the angle of the distance vector r^, 
with rji = rji (cos aji, sin a ji). From f)2.5p . one finds the two different spatial regions 
of approaching and receding particles, respectively, where the relative velocity has a 
different sign for fixed ipi and ipj in the interval (fi < ipj < ipi + 2ti: 



Vji > for fi_j£t < a „ < <£i_jh _|_ n m0 ving away 
Vji < for _|_ vr < a .. < n±£i _|_ 2-ir approaching. 



Hence, a particle located in the half-sphere in clockwise direction from the mean angle 
(tpj + <fi)/2 approaches the focal particle i, whereas particles located in the other half- 
sphere anticlockwise from the mean angle are moving away (see figured]). Please note 
that for ifj = (fi the social force vanishes as Vji = according to (12. 5p . 

Four different regions in the (/i m , /x a )-parameter space are distinguished [28] : 

• Pure Repulsion: /i m < 0, \i a < 0; 

• Escape/Pursuit: /i m > 0, \i a < 0, i.e. attraction to particles moving away, repulsion 
from particles coming closer; 
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Figure 1. Visualization of the spatial regions for approaching and receding self- 
propelled particles for a binary interaction of the ith particle with velocity vector v, 
(heading angle </?,), and a neighbouring particle j with velocity Vj (heading angle ifj) 
within its sensory range l s . The center of the circle corresponds to the position of the 
focal particle i. The dotted line, determined by the mean angle ((pi + ifj)/2, represents 
the border between the two distinct spatial regions (half-discs) corresponding to 
approaching and moving away of individual j. If the jth particle is located above 
the dotted line in the red region, the two particles are coming closer (approaching). 
Otherwise, if the neighbouring particle is located in the blue region, the two particles 
move away from each other (as shown in this example). 



• Head on Head: /i m < 0, /i a > 0, i.e. attraction to particles coming closer, repulsion 
from particles moving away; 

• Pure Attraction: /i m > 0, fi a > 0. 

Some typical spatial snapshots obtained from individual-based simulations of (12. ip are 
shown in figure [2j 

The above microscopic model corresponds to the ones studied in [271128]. if 
the distance dependent interaction strengths fi r (rji) and n a ,miXji) are assumed to be 
constant and the forces are rescaled by the number of particles within the interaction 
area of the focal particle. 

3. Kinetic Description 

In this section, we derive a kinetic description for the above individual-based model 
(12. ID . For this purpose, we introduce the iV-particle probability density function (PDF) 



which determines the probability to find particle i at position T*j moving into the direction 
(fi (i = 1,2,..., N) at time t. It is normalized with respect to integration over all 
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Figure 2. Spatial snapshots of the microscopic model for different interaction 
strengths after the system relaxes towards a steady-state: (A) disordered, homogeneous 
state for [i a = \x m = —0.6 (pure repulsion), (B) diffuse collectively moving bands for 
fi a = —0.6, fi m = (only escape), (C) collectively moving bands for /i a = —0.6, 
fx m = 0.6 (escape and pursuit), (D) dense, collectively moving cluster for fj, m = 0.6, 
fx a = (only pursuit), (E) cluster state for (j, m = fi a = 0.6 (pure attraction), (F) 
disordered, homogeneous state for [i m = —0.6, fi a = 0.6 (head on head). The velocity 
vectors of individual particles are coloured corresponding to their polar direction of 
motion according to the inset of (A). The inset in (D) indicates the positions of the 
snapshots in the (/i a , /U m )-parameter plane. Other parameters: N — 4000, L = 40, 
D v = 0.02, l s =l,l c = 0.2, fi r = 5, s = 0.25. 



positions and angles 

N , . 



f[(^J <%J Pjv(n, <pi) r 2 , <P2',.-.', r N , <p N ; t) = 1 . (3.2) 

In agreement with ( 12. ip . we can write down the Fokker-Planck equation (FPE) for the 
dynamics of the PDF as follows. 



dP N 



Self-propelled particles with selective attraction-repulsion interaction 



7 



From the linear FPE above one can derive an evolution equation for the marginal PDF 

dipj ) P N (r u <pi, 7- 2 , ip 2 \ ■ ■ ■ ; r N , <p N ; t) (3.4) 

by integrating (13. 3p over the degrees of freedom of particles j ^ i. In what follows, we 
assume that correlations between particles can be neglected. Therefore, the iV-particle 
PDF shall factorize, i.e. 

N 

Pn{?\, tpx\ r 2 , tf>2\...\ r N , Lf N ; t) = JJ P(r h <p it t). (3.5) 

i=i 

By this means, one obtains an effective one-particle description 

dP(ri,<pi,t) N-l d D V &P 

— = -s V ri {e i)h P) — (F v P) + -^ ri r^; » 3 - 6 

dt s d(fi * s^ dtpf 

where the force F v is given by the following integral over P(ri, (pi,t): 

"ifi+2-K 



F v (r h <pi,t) =2s f dipj sin ( — ^ J > ) / d r ji r ji 

J(fi V ^ / Jlc 



Vm(rji) I ^ daji sm(aji - </?i) sin ( ^ - — - ) P(r, + r^, pj, t) 



2 



-A*o(»ji) da i* sin(<x,-j - sin ^ — ^ ^ J - ay^ P(r< + r^, <pj, t) 

pifi+2-K nl c p2ir 

-I dtpj dr^rji da j{ /^(r^) sin^ - (pi)P(ri + r jit <pj,t) . (3.7) 
JiPi Jo Jo 

The factor N — 1 in (13.61) arises, because (I3.3P was integrated over the positions and 

angles of N — 1 identical particles. In other words, the focal particle can interact with 

N — 1 neighbours. For the same reason, the particle index % is omitted henceforward. 

Please note, that a nonlinear FPE (13. 6p for the dynamics of P(r,ip,t) was derived 
from a linear FPE (13. 3 p governing the dynamics of P^. Equation (13. 6p is nonlinear 
since the interaction force F v depends on the probability density P(r,ip,t). In this 
sense, assumption (13. 5p is a mean-field approximation: a single particle is affected by a 
force due to its own PDF (13.71) . 

By introducing the one-particle density function p(r,tp,t) = NP(r,tp,t) we can 
eliminate the factor (N — 1) « N from (13. 6p . Accordingly, p(r,tp,t) is interpreted as 
particle density obeying the nonlinear Fokker-Planck equation 

Mr,<P,t) , , d F ip {r,<p,t)p d 2 p 

dt = ~ SoVr {ehp) ~ dip 7 + SJW ' ( } 

One solution to the equation is the uniform distribution (spatially homogeneous, 
disordered state), i.e. 

P (0) = £ , (3-9) 

where po — N / L 2 is the homogeneous particle density and L denotes the linear spatial 
dimension of the two-dimensional square-shaped system. 



Self-propelled particles with selective attraction-repulsion interaction 



8 



3.1. Stability Analysis in Fourier Space 

In order to analyze the Fokker-Planck equation ( 13. 81) analytically and especially to solve 
the integral (13. 7p . it is convenient to work in Fourier space with respect to the spatial 
coordinates r and the angular variable tp. The dynamics of the Fourier coefficients 

g n (k,t) = J d 2 r ^ d<p P (r,<p,t)e ikr+i ^ (3.10) 

reads 

d 9n{k,t) tk x s kyS r,2 D <PZ 
7^ 2 — K9n+1 + 9n-l) — {9n+l ~ 9n-l) ~ ~n ~^9n 

oo „ 

Yl I d 2 q gj (q,t)k n ^(k-q,t), (3.11) 

;'=— oo 

where the Fourier coefficients of the force (13 .7p are defined as 

K n (k,t)= [ d 2 r d<pFJr,<p,t)e ikr+in *. (3.12) 



(27r) 3 s 

j=— oo 



Obviously, the Fokker-Planck equation turns into an infinite hierarchy of equations 
(13. lip in Fourier space. 

In this section it is shown, how the stability of the uniform distribution (13. 9p can be 
analyzed in Fourier space. Furthermore, some difficulties associated with the derivation 
of hydrodynamic equations, as done in section HJ are outlined. A similar approach was 
used by Chou et al in order to achieve a kinetic theory for self-propelled particles with 
metric- free interactions [2U| . 

For the stability analysis, the Fourier coefficients K n (k,t) of the force (13. 7p are 
required. Due to the fact, that the force F v (r, <p, t) depends on the integral over p(r, <p, t) 
(13. 7p . the Fourier coefficients K n (k, t) of the force can be written as a linear combination 
of the Fourier coefficients g r (k,t) 

oo 

K n (k,t)= K 1hr (k)g r (k,t) (3.13) 

r=— oo 

using an infinite dimensional matrix K n ^[k). More details on the calculation of the 



matrix elements are given in Appendix A Consequently, the nonlinear Fokker-Planck 
equations reads 



dg n {k,t) ik x So k y SQ . 2-^<p 
-Q t y~ ( 9n+1 + 9n -^ 2~ ^ 9n+l ~ 9n ~ 1 ' = ~ U ~sj- 



oo oo „ 

+ 7^^ H I d 2 qg 1 (q,t)g r (k-q,t)k n _,, r (k-q). (3.14) 

^ ' " j=— oo r=— oo 

Let §n\k,t) be a solution of the equation above. The dynamics of a small 
perturbation Sg n (k,t) = g n (k,t) — cjn\k,t) is determined by the linearized equation 

d5g n (k,t) ik x s k y s in 
{dg n+1 + og„-x) — {og n+ i - og n ^i) = -n -^-og n 



at 2 2 Sq (27r) ci s 
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oo oo „ 

J2 E / ^ (%(qM°\k-q,t)+gf\q,t)6g r (k-q,t))k n ^ r (k 

j=—oo r=— oo 

The Fourier transform of the uniform distribution (13.91) reads 



(3.15) 



g^(k)^(2nfp^5(k)S n>0 . (3.16) 

It is isotropic, i.e. it is invariant under rotations, and constant in time. Therefore, 
one can choose without loss of generality k = (k,0). Inserting ( 13 . 1 6[) into (13. 15[) and 
using that the force F v vanishes in the homogeneous, disordered state, yields an infinite 
dimensional system of linear differential equations: 

dSg n (k,t) _ ^ {iksQ 2 D V mp (0) ~ 



dt 



{ — - (S n +i,r + S n -i,r) - n 2 —!f5 n>r H K n>r (k) > Sg r (k,t) . (3.17) 

I A Sr. SO J 



Hence, the stability of the spatially homogeneous, disordered state is determined by the 
eigenvalues of the matrix 

ikso inp^ ~ 
M n:r = — — (d n+ l, r + f)„-l,r.) - u — 5 n , r H . (3.18) 

In order to calculate the eigenvalues numerically, it is necessary to find an appropriate 
closure of the infinite dimensional system of linear equations. Here we assume, that a 
critical h > exists, such that 5g n = for \n\ > n holds (cf. [29]). That assumption is 
reasonable, because the nth Fourier coefficient is strongly damped, since 

dS ^ t] oc-n^Sg n (k,t). (3.19) 

L/t/ Sq 

Using this approximation, the stability of the spatial homogeneous state can be studied. 
Please note, that only two approximations were used: First, correlations between 
particles are neglected (13. 5p (mean- field) ; Second, the truncation of the hierarchy of 
equations (13.171) . In principle, it is possible to analyze the stability for arbitrary 
wavenumbers k = \k\, notably, there is no restriction to small k. Figure |3] shows the 15 
largest eigenvalues for different parameters as a function of the wavenumber. 

Unfortunately, the analysis discussed above does not supply any analytical criteria 
for the stability of the solution p(°\ such as critical microscopic model parameters. 
Those will be derived in section |4] from the hydrodynamic theory. Besides, p^ is the 
only solution to the Fokker-Planck equation known analytically, which stability one can 
investigate. 

In this context, we would like to emphasize one difficulty which occurs if the focus 
is put on the dynamics of the lowest Fourier coefficients on large length scales, i.e. small 
wavenumbers k: The matrix elements K n ^ r (k) involve Bessel functions of the first kind 



of k (see Appendix A and [29J) which are oscillating functions with alternating Taylor 



coefficients [30]. The approximation 



^ (-l) m /X\2™+" 



m=0 
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Figure 3. Eigenvalues of the truncated matrix M„. r f|3 . 1 8[) as a function of the 
wavenumber k. The 15 largest eigenvalues are shown, i.c n = 7 Fourier coefficients 
are considered. Different scenarios are distiguished: (a) fi m = —2, /i a = — f: the 
spatial homogeneous solution p 1 - ^ is stable; (b) \x m = 5, fi a = 5: long-wavelength 
instability; (c) fi m — 2, fi a — — 1: homogeneous state (k = 0) becomes unstable first; 
(d) ii m = —3, fi a = —3: destabilization of a finite range < k m i n < k < k max of 
wavenumbers. Other parameters: D v = 0.1, po = 1, fi r = 0, l c = 0, l s = 1. 



will only be valid in the vicinity of x — 0. For large x, J Vj n( x ) tends to either plus 
or minus infinity, but never to zero. This issue will be important for the derivation 
of hydrodynamic equations in the next section and is restricting their validity to small 
wavenumbers and large length scales, respectively. 

4. Derivation of Hydrodynamic Equations 

In this section, a hydrodynamic description of the many-particle system is derived 
directly from the one-particle FPE (13. 8ft . For this purpose, it is convenient to work 
in Fourier space with respect to the angular variable tp. Both, the particle density 
p(r, f, t) and the force F v (r, tp, t) (13. 7p are 2-7r-periodic functions in tp and are therefore 
expanded in a Fourier series as follows: 

1 °° 

p{r,tp,t) =— Ur,t)e- in \ (4.1a) 

n=—oo 
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.. oo 

Fy(r,<p 7 t) = — Qn(r,t)e- in <?. (4.1b) 



2tt 

n=— oo 



The Fokker-Planck equation in Fourier space is obtained by multiplying (13. 8p with e mLp 
and integration over ip e [0, 2ir]: 

% + (Yfn-i + Y*/ n+1 ) =~T, fi - n *^rf» ■ ( 4 - 2 ) 

In (14. 2p . the complex derivative V = 0.5 (d x + id y ) was introduced. The Fourier 
coefficients 

f n (r,t)= d<pp(r,<p,t)e in * = p(r,t)(e* n *) (4.3) 
Jo 

obey the symmetry f*(r, t) = f~ n (r, t), since p(r, <p, t) is real- valued. The lowest Fourier 
coefficients are related to macroscopic physical quantities, namely the marginal particle 
density p(r,t) and also the momentum field w(r,t) = (w x ,w y ) via 

p(r,t) =f {r,t), (4.4a) 

to«(r ) *) = Y(/i(r,*) + /-i(r,*)) , (4.46) 

"v(r,0 = ^(A(r,0-/-i(r,0) ■ 
The symmetric temperature tensor 

f = ( (cosf 2 ) - (cosip) 2 (sin cp cos ip) - (sm(p) (cos 99) \ g 

I (sin <^ cos — (sin</j) (cos<^) (sin(/? 2 ) — (siny?) J 

as defined in [23l[2ni[2B] , is basically related to the second Fourier coefficients: 

1+ A(r,Q+^(r,«A_/ , / 1 (r,«)+/ 1 (r, t )y 
2p(r,t) ) \ 2p(r,f) ^ 

= A(rt)-/_,(r,«) _ jkr.tj)' - (Ur.t)) 1 
" 4 V (r,() 4i(p(r,i)) 2 

T " = 2( j + ( WJ) J ' (4 ' 6C) 

In the following, the behaviour of these quantities is deduced from (14.21) . Therefore, 
the Fourier coefficients of the force Q n (r,t) are approximately calculated by expanding 
p(r + rji,ip,t) into a multidimensional Taylor series for small |r\,j| [31], substituting 
(14.1 aj) into ( 13.71) and evaluating the remaining integral. Here, we consider all terms up 
to the third order of the Taylor expansion. This approximation holds for spatially slowly 
varying densities, i.e. large system size compared to the interaction radius l s . 

Again, the Fokker-Planck equation turns into an infinite hierarchy of equations 
(14. 2 p in Fourier space. An appropriate closure scheme is used in order to consider only 
the first Fourier coefficients. Close to the order-disorder transition, the mean speed of 
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the system is small, i.e. \w\ /(sop) ^ e > where e is a small number. Furthermore, the 
dynamics of the Fourier coefficients (14.21) suggests, that fi is larger than /*. with \k\ > 1, 
because of the damping term in (14. 2ft proportional to k 2 , leading to the scaling relations 

Voce, <9 t oce, f k oc . (4.7) 

In [20] it is argued, that the scaling ansatz of the temporal and spatial derivatives reflects 
the propagating nature of the system. This closure scheme described above was already 
used in [T71[20l|22l[3T] and in [21] in an analogous way for nematic particles. A system 
of three nonlinear partial differential equations is obtained by keeping all terms up to 
the order e 3 . 
dp 



dt 

dfi 
dt 



s YA* + V7i 



32 



(6p - 6) h - £1/2/1* + 6 - 7T& hTp + t^p - so Th 



9 



32 



45 



32 



+ -w&P ~ &P -so VP- tt^/iY/i 



9 



64 



45 



2 



p a/i + 2vvr - 2/r Y 2 P + (e 7 - y*(y>) 



5/2 

9* 



A + &/? - + %p] y/i + f f e 3 - & ) A yp - ^py 2 



2 45" r y 



The following coefficients are introduced: 



(4.8 a) 



(4.86) 
(4.8c) 



£1 = 


7r r a 

- / drjiTji [/i m (r>) - Pa{rji)\ , 

J l c 


(4.9a) 




— I dr ji r2 jiPr(rji) , 
s o Jo 


(4.96) 


£3 = 


1 f ls 

- / dr i4 rj [/x ro (rj-i) + Pa(?>)] , 


(4.9c) 


£4 = 


8 J dTji ^ ~ ' 

1 > 
s\ 


(4.9 d) 


£5 = 


(4.9 e) 


£e = 


IT f lc 

2s~olo drjir ^ ifir ^ ' 


(4-9/) 


6 = 


16 /"' s 

— / d rji Tji WiTji) + Mafai)] • 

9tt 


(4.9<?) 



Equation (14.8 aft is simply the continuity equation representing the conservation of the 
particle number N. 
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4-1- Hydrodynamic Limit 



Another simplification of the system f]4.8)) is reached by adiabatically eliminating /2, i.e. 
dtf2 ~ 0, where 



s , 64, 



45 



32 



9 



(4.10) 



follows. Please note, that £5 oc > was assumed. By identifying Fourier coefficients 
with physical quantities (14.4j) . familiar hydrodynamic equations are obtained. 
dp 



-V ■ w 



(4.11a) 



at 

-^r = (Ai - 771 |™| 2 ) to + A 2 Vp + A 3 V(Ap) + X A Aw + A 5 V • (V • w) + rj 2 (w - V) • w 



+ V3 



1 2 

-V |tw| — w ■ (V • it/ 



- r] 4 (w ■ V) • Vp + 775™ ■ Ap + 77 6 tt; ■ ((to • V)p) 

+ 7l w ■ (Vp) 2 + 72 (V ■ id) ■ Vp + 73A, • Vp + 74M™ • Vp 
+ 75 [2(Vp-V)-Vp-(Ap)-Vp] 
Here, the matrices 



(4.116) 



d x w x 



\ (dyW X + d X Wy) dyWy 







d X Wy - dyW X 



8yW X ~ 8 X Wy 



(4.12 a) 



(4.126) 



were introduced. The transport coefficients {Xi,^,^} fIB.lj) as functions of & (14. 9 p are 
listed in Appendix B The hydrodynamic equations (14. lip are comparable in structure 
to the theory of Toner and Tu [i~3|[T4"] but contains additional gradient terms comparable 
to [IS]. 

Let us begin the analysis by neglecting the spatial derivatives in (14. lip and analyzing 
the fixed points of 

— = (Ax-r^H 2 )™. (4.13) 

Apparently, two spatially homogeneous solutions exist: First, w = corresponding to 
a disordered phase and vanishing center of mass velocity. Second, 



w 1 = \ —e. 



(4.14) 



where e denotes an arbitrary unit vector reflecting the isotropy of the system. For 
Ai < 0, only the fixed point w exists and is stable against spatially homogeneous 
perturbations. For Ai > 0, Wi is a second fixed point corresponding to an ordered phase 
with nonzero mean-speed (swarming phase). Please note, that 771 is always positive. In 
Ai = 0, i.e. 

4LL 



d.Tj % Tj i \fJ"rn yj i) l^aiXjif) 



vrp sg ' 



(4.15) 
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the system undergoes a supercritical pitchfork bifurcation. By introducing the effective 
coupling strength 



l> 



drjirji(^fi m (^rji^j /i a (rjj)) 



lc 



together with the critical coupling parameter 
4A, 



He 



vrpoSo 



(4.16) 



(4.17) 



one can rewrite (14. 14)) close to the critical point < /i — p c <C p c in the leading order 
as follows 

;9 



w l 



soPo 



V - Pc 

Pc 



(4.18) 



where one can read off the mean-field exponent /3 = 1/2. Hence, the order- disorder 
transition is continuous in the mean-field analysis. 

Equation (14.151) defines the transition line in the (/x m ,/x a )-parameter space. 
Interestingly, only the integrated interaction-strengths are important, the exact 
functional dependence on the distance between two particles p m ^ a = fJ>m ta (rji) does not 
matter. Furthermore, collective motion is possible in both situations, pure attraction to 
particles moving away (pursuit: fi m > 0, \i a = 0) and pure repulsion from approaching 
particles (escape: p m = 0, p a < 0), respectively. 

Up to now, only the spatial homogeneous system was analyzed. In the following, the 
stability of the disordered, homogeneous solution wq = against spatial perturbations 
is considered. In order to investigate the stability of w Q , the ansatz 

p(r,t) =p + 5p{r,t), (4.19a) 
w(r,t) = 6w(r,t) (4.196) 

is inserted in (14. lip and the resulting equations are linearized in the perturbations which 
are assumed to be small. The linearized equations read 
85 p 



dt 

d5w 



-V -Sw. 



(4.20 a) 



XiSw + X 2 VSp + A 3 V(A5p) + XiASw + A 5 V ■ (V ■ 5w) . (4.206) 



Inserting the exponential ansatz 

5p(r , t) = ( p exp (at + ikr) , 
8w(r, t) = C, w exp (at + ikr) 

yields the growth rate a depending on the wavenumber k 

Ai — X^k^ , 

Ax - (A 4 + A 5 ) k 2 ± V / [Ai-(A 4 + A 5 )^] 2 



Ifcl 



1 

*2,3 = ^ 



Ak 2 (A s 



A, k 2 ) 



(4.21a) 
(4.216) 

(4.22 a) 
(4.226) 



Provided all eigenvalues are less than zero, the disordered spatially homogeneous state 
is stable. The first eigenvalue a\ is equal to Ai for k = 0, i.e. the stability of Wq is 
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determined by the sign of Ai. That is the destabilisation of the homogeneous state 
as discussed before. For A4 > 0, the first eigenvalue o~\ is monotonically decreasing. 
However, the eigenvalue o\ is monotonically increasing for A4 < 0. In this case, the 
hydrodynamic equations (14. lip loose their validity. This unphysical behaviour is due to 
the restriction to third derivatives when the integral ( 13.71) was calculated. Furthermore, 
it is problematically that the Taylor coefficients have alternating signs so that Taylor 
polynomials converge slowly as argued at the end of section |3~T1 To predict the behaviour 
of the system in that case, one needs to consider higher order terms. Unfortunately, 
the number of terms that has to be considered to be consistent with the scaling ansatz 
(14.7j) . grows rapidly. 

Nevertheless, it is possible to predict the stability of w by analyzing those 
eigenvalues, which tend to zero for small wavenumbers (hydrodynamic modes |14j). 
For small wavenumbers, the eigenvalue 0*2 is expanded in a Taylor series: 

a 2 ~ k 2 + Xl + A ' As ~ y 2 (A4 + As) k 4 + O(k') . (4.23) 
Ai X 1 

Suppose, w is stable against spatially homogeneous pertubations, i.e. Ai < 0. In 

that case, w will become unstable, if A2 > 0. Considering (14.111) . this instability is 

reasonable: for A2 > 0, inhomogeneities in the particle density will lead to a flow in the 

direction of the density gradients while for A2 < density inhomogeneities decay due 

to inverse flows. The amplification of density gradients will lead to an agglomeration of 

particles (cf. figure [2j5) . This hypothesis is supported by numerical simulations of the 

microscopic model, as shown in the next section and in [28]. A similar instability due 

to the gradient term in ( 14. lip is also known from other active matter systems [221132] . 

The critical line is given by A2 = 0, i.e. 

f i ^rjWrjil+^y^^ + ^l dr^rj^fai)) . (4.24) 

The stability analysis of the ordered state Wi using the full hydrodynamic equations 
is much more demanding. It is shown in [IT] that a long wavelength instability can occur 
in the ordered regime. Furthermore, giant number fluctuations are expected to emerge 
because of the structure of the equations (14. lip [i~5l[T6] . 

Finally, the prediction of the hydrodynamic theory derived in this section is 
compared to the predictions of the stability analysis using the kinetic description in 
section 13.11 Figure H] shows the eigenvalue spectrum by the two theories, whereby 
the hydrodynamic description ( 14. lip can be understood as an approximation to the 
full system (13.111) . The hydrodynamic theory describes the dynamics of the many- 
particle system on large length scales and small wavenumbers, respectively. Thus, the 
eigenvalues coincide in the limit k — > 0. For large k, deviations occur due to the 
approximations made, namely: 

• restriction to third derivatives; 

• scaling ansatz and closure scheme (14 .7p : 

• restriction to the lowest Fourier coefficients. 
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Figure 4. Eigenvalues as predicted by the hydrodynamic theory (|4.11j) (red lines) 
and the kinetic description p. lip (black crosses) which determine the stability of 
the disordered, spatially homogeneous state: = pq/{2-k)\ p(r,t) = po = N/L 2 ; 
w = w = 0. For small wavenumbers k — > 0, the eigenvalues coincide, i.e. the 
hydrodynamic theory is valid on large length scales. Parameters: (a) p m = —2, 

p a = -1; (b): p m = 2, p a = ~1; D v = 0.1, p = 1, p r = 0, l c = 0, l s = 1. 



Nevertheless, the hydrodynamic theory yields important insigths on the macroscopic 
behaviour of the system. The symmetry breaking, i.e. the existence of an collective 
motion mode is predicted. Moreover, the stability of the homogeneous solutions to the 
hydrodynamic equations can be analyzed both analytically and numerically. Finally, it 
allows to obtain analytical results on the order-disorder transition line as well as the 
occurrence of clustering. 



5. Comparison with Numerical Simulations 



In order to analyze the stability of the disordered, homogeneous state, we have performed 
systematic numerical simulations of the individual-based model (12. ip . The degree of 
collective motion was measured using the time averaged polar order parameter 



1 



N 



(5.1) 



Here (-)t represents a temporal average. ($) f = 1 corresponds to perfect orientational 
order with all agents moving in the same direction, whereas a vanishing corresponds 
to a completely disordered system. Please note that = can only be observed in 
the thermodynamic limit (N — > oo). In a finite, disordered system we will measure a 
small, but finite (3>) t <; due to finite size fluctuations of the order 1/y/N. 

In order to measure the deviations from a homogeneous state, we have subdivided 
the simulation domain into square cells of size l s x l s set by the sensory range. We used 
this spatial subdivision to calculate the spatial entropy function 



S 



Pi l °SPj = - 



— log — , 

N S N' 



(5.2) 



Self-propelled particles with selective attraction-repulsion interaction 



17 



where the summation occurs over all occupied cells of the grid with rij > 1 (nj is the 
number of particles in the jths cell). This allows us to define the following spatial order 
parameter 



with S max being the maximal value of the spatial entropy corresponding to a 
homogeneous distribution of particles. ($) t = corresponds to a perfectly disordered 
state, whereas > indicates a spatially inhomogeneous distribution of particles 
(clusters, bands), where ($) t = 1 corresponds to the extreme situation where all particles 
are located in a single cell. 

In order to test the stability of the disordered state, we have averaged the two 
order parameters over a time interval At = 1000 after an initial time = 1000. 
Please note, that for certain parameter values it is possible that the system has not 
reached a stationary state after t = 1000. However, this initial time is sufficient to 
account for deviations from the homogeneous disordered state, which we are interested 
in. Accordingly, we are not interested in the actual stationary values of and (&) t . 
Thus, we set the upper limit of the colour bar for both order parameters in figures [5] and 
[6]to 0.2. In consequence, all regions of parameter space, where (<E>) f > 0.2 or (\&) t > 0.2 
will appear white in figures [5] & [61 

The stochastic differential equations of the microscopic model were integrated with 
periodic boundary conditions using the stochastic version of the Euler algorithm with 
a numerical time step dt = 0.01. The initial condition for all simulations was the 
disordered, spatially homogeneous state. 

For simplicity, we consider only constant interaction strengths (independent of the 
distance): n m ,aijji) = /i miQ = const. Here, we focus on the analysis of the stability 
of the disordered, homogeneous solution with respect to variations of the interaction 
parameters fi a and fi m . The analytical predictions of the hydrodynamic theory for 
the onset of instability of the disordered solution are represented by two intersecting 
critical lines in the (/i a , /i m )-plane perpendicular to each other (figure [5]). The first one 
(/j, a ~ fi m ) corresponds to the orientational instability and onset of collective motion, 
cf. (14.151) and line (1) in figure whereas the second one (fi a ~ — /i m ) corresponds to 
the density instability associated with structure formation due to effective attraction 
between particles, cf. (14.241) and line (2) in figure El The homogeneous, disordered 
solution is predicted to be linearly stable only "below" both critical lines. 

The emergence of collective motion as well as the destabilization of the homogeneous 
density distribution is mostly confirmed by the numerical individual-based simulations. 
In specific parameter regions, deviations from the hydrodynamic theory are observed. 

In particular, we observe a high degree of collective motion in the escape and pursuit 
region and strong deviations from the homogeneous spatial distribution of particles in 
the pure attraction regime as well as in the escape-pursuit regime (in particular in the 
pursuit dominated regime for \fj, m \ > |/i a |). This also confirms our previous results of 
comprehensive numerical investigations of related individual-based models [281133]. 




(5.3) 
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Ha Ma 

Figure 5. Comparison of numerical simulations of the individual-based model with 
the predictions of the hydrodynamic theory on the stability of the disordered solution 
in the (fx m , /i a )-plane for different densities: po = 10.0 (L = 20, top) and po = 2.5 
(L = 40, bottom). Left: Average orientational order parameter ($)*; Right: Average 
spatial inhomogeneity order parameter (^)t- The solid blue (gray) lines show the 
critical lines for the instability of the disordered solution. Line (1) corresponding to the 
orientation instability (|4.15[) . whereas line (2) corresponds to the clustering instability 
621. Other parameters: D v = 0.06, N = 4000, fj, r = 5, l c = 0.2, ; s = 1, s = 0.25, 
dt = 0.01. 



Disagreements between simulations and the theoretical prediction (non-vanishing 
orientational order and/or clustering in simulations below the two critical lines) appear 
predominantly close to the intersection of the two critical lines and are stronger at low 
densities. They might be associated with the mean-field assumptions used in order to 
derive the hydrodynamic theory. On the one hand, at low densities the assumption of a 
continuous density of neighbours is strongly violated. On the other hand, correlations 
between particles may play an important role in the respective parameter range (likewise 
for high densities). Thus, the factorization of the N-particle PDF into a product of one- 
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Figure 6. Results of numerical simulations on the stability of the disordered solution 
in the (/x m , /z a )-plane for vanishing short-ranged repulsion (l c = 0). The spatial 
order parameter (\&). shows emergence of structures below the critical line (2). Other 
parameters: N = 4000, l s = 1, s = 0.25, D v = 0.02, L = 40. 

particle PDFs f 13 . 5 p leads to a questionable approximation in that case. 

Another possible explanation for disagreement between theory and simulation is a 
breakdown of the homogeneous, disordered solution due to finite amplitude instabilities 
at parameters where this solution is still linearly stable. 

For weak (or vanishing) short-ranged repulsion (l s ^> l c or ii r <C |j«m,a|)> we observe 
inhomogeneous states without polar order far in the head-on-head regime (/x m < 0, 
fx a > 0), clearly below the critical line predicted by the hydrodynamic theory as shown 
in figure [6] This instability was missed in the previous study of the model [28]. We were 
able to confirm it, using our kinetic approach (see section [3]) by positive eigenvalues 
of the matrix determining the stability of the disordered solution at the respective 
parameter values. A close inspection of the dynamics of the individual-based model in 
this parameter region reveals the emergence of dense nematic filaments with particles 
moving in an anti-parallel fashion within the filament, whereby approximately half of 
the particles moves in either direction along the filament (see figure [7]). 

The hydrodynamic theory derived in section H] is not able to account for this kind of 
structures as it considers only the (polar) momentum field and not the nematic director 
field as a coarse-grained variable. An extension of the hydrodynamic theory by including 
an equation for the nematic director field may allow for identifying instabilities due to 
the onset of nematic order, however it goes beyond the scope of this work. We refer the 
reader to some recent theoretical works on active nematics [2T |l34ti37] . 
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Figure 7. A sequence of snapshots of one simulation for different times visualizing the 
formation of nematic filaments for fi a = 0.6, /i m = —0.4. Other parameters: N = 4000, 
l s = 1, s = 0.25, D v = 0.02, L = 36. 



6. Discussion 



In this work, we have derived and analyzed a kinetic and hydrodynamic description of 
self-propelled particles with selective attraction-repulsion interaction. 

At first, we derived a kinetic description of the system based on the Fourier 
transform of the probability density function. The corresponding system of equations 
for the successive Fourier modes can be used for efficient numerical analysis of the linear 
stability of solutions of the nonlinear Fokker-Planck equation. We have analyzed the 
stability of the disordered, spatially homogeneous solution. In addition, we have shown 
that the integration over the social forces yields Bessel functions of the first kind, which 
enter the matrix elements of the corresponding linearized system of differential equations 
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in Fourier space. Due to the alternating Taylor coefficients of the corresponding 
Bessel functions, a closure approximation, corresponding to a finite order expansion, 
immediately leads to unphysical divergences at large wavenumbers k. 

Furthermore, we have derived a hydrodynamic theory by truncating the "small- 
wavenumber-expansion" at the third order. The resulting hydrodynamic equations 
are in agreement with the generic Toner and Tu theory of active matter but contains 
additional gradient terms due to the selective attraction-repulsion interaction. Our work 
establishes a direct link between the microscopic parameters of the individual-based 
model and the macroscopic parameters governing the behaviour of the coarse-grained 
hydrodynamic variables (density and momentum fields). A comparison between the 
eigenvalues obtained from hydrodynamic theory and the kinetic description, which takes 
higher orders into account, allows us to assess the validity of the hydrodynamic theory. 

We performed extensive simulations of the microscopic model based on stochastic 
differential equations focusing on the (/x m , /i a )-plane in the parameter space. The 
critical lines obtained from the hydrodynamic theory, where the disordered, spatially 
homogeneous solution becomes unstable, show good agreement with the numerical 
results at sufficiently high densities. However, at certain parameters clear deviations 
appear, as for example in the vicinity of the crossing of the two critical lines at low 
densities. 

This leads us to the important question of the validity of the approximations 
made in order to derive the coarse-grained theory. One common assumption is that 
of molecular chaos, which allows to factorize the N-particle PDF into a product of 
N one-particle PDFs. In models with collision-like interactions, the approximation is 
supposed to work best at low densities, where the mean-free path of particles between 
interactions is large [TTl[T9] . However, in our case, we observe large deviations at low 
densities. This contradicting effect may be due to an approximation required to evaluate 
the integrals over the social forces. It relies on a Taylor expansion of the one-particle 
density function around the position of a focal individual. Thus p(r,<p,t) is assumed 
to be continuous and differentiable, an assumption which is likely to be violated at 
low densities. Similar arguments have been put forward also in (22]. The impact 
of individual (angular) noise on the mean-field assumption may also be not straight 
forward. Intuitively, one would argue that uncorrelated individual noise terms always 
decrease correlations between interacting particles. However, for self-propelled particles, 
angular noise leads to a stronger localization of particles [SEIEH], thus in principle also 
to a prolonged interaction between neighbours, which may in principle enhance multi- 
particle correlations. 

Furthermore, the systematic comparison of the numerical results with the analytical 
prediction revealed an unexpected additional instability. It corresponds to the 
emergence of filamentous structures with nematic order. Onset of nematic order is linked 
to the dynamics of the second Fourier amplitude, which was adiabatically eliminated 
in order to derive the hydrodynamic theory. Thus, the hydrodynamic theory cannot 
account for this instability, but we were able to confirm it by numerical evaluation of 
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the linearized kinetic equations derived in section El 

So far, hydrodynamic equations of active matter were derived directly from minimal 
microscopic models of self-propelled particles with velocity-alignment. Here, we show 
that it is also possible to derive such equations for a more complex model of self- 
propelled particles with selective attraction-repulsion interaction and establish a direct 
link between the microscopic and macroscopic level of description. The model exhibits 
a large variety of different phases. Thus we believe, it might be not only of interest 
from the biological point of view, as an alternative to models including explicit 
alignment of individual agents, but offers also interesting playground for the study 
of self-organization, pattern formation and phase transitions at far-from-equilibrium 
conditions. 

Appendix A. Calculation of K p>n (13.131) 

In this section, it is sketched, how the matrix elements K PtTl (k), defined by 

// 1 2-7T 
d 2 r J d<pF(r,<p,t)e ikr+ipip , (A. la) 

oo 

K p {k,t)= K P , r (k)g r (k,t) (A. 16) 

r=— oo 

and used for the stability analysis of the spatial homogeneous, disordered solution ( 13.91 ) 
of the Fokker-Planck equation (I3.8P in section 13.11 are calculated. K pn is obtained by 
inserting the inverse Fourier transform of the one-particle PDF 

oo - 

P(r , <P,t) = j-~ 3 I d ' k 9n(k, t)e~* kr - m * (A.2) 

^ ' n=— oo 

into (13.71) and solving the remaining integral: 

1 OO „ f / \ 

F^r^t) = 7^3 E / d 2 kg p (k,t)e' ikr J d^e"*" j^sm [^^j 



p=—oo 



f , , , [~ + \, - lkr ■ , ^ ■ (<P + <Pj \ 

J drjiTji n m {rji) J v+v . dae J * sm(a - ip) sin ^ — 2" ~ y 

~Va(rji) J v+ l. da e~ ikr ^ sin(a - <p) sin { ^^ 3 - aj 

<!/>•/>• //,.(/>•) j da sin(a - ^)e~ ikr ^ j . (A.3) 

In order to perform the integration, the exponential function e~ lkVji is expanded into its 
Taylor series 

OO / .X g 

e~ ikr " =Y, - Z T~ cos ° (x-a), (A.4) 



=0 s! 
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where k = k (cosx, sinx) and rji = (cos a, sin a). Due to the isotropy of the uniform 
distribution (I3.9p . which stability is under consideration, x = is chosen without loss 
of generality. Furthermore, it is convinient to use the following identity for s6N [30] : 



cos (a 



) = h E (f) cos (( s - 2/ )«) • 

1=0 \ ' 



(A.5) 



Performing the integral yields the force expanded into a Fourier series and as a linear 
combination (lA.lfc|) of the Fourier coefficients g r (k,t), where one can read off the 



elements of the infinite dimensional matrix K Pjn (k): 

K Pt n(k) = —^r- [(p m ,o{k) - Hafl(k)) (5p,l5 n ,l ~ 8p ,-l5 n -l) 

+ {p m ,2(k) - fJ>a,2{k)) (5p-l5 n ,l - S Pt i6 n -i + (5 Pt2 - S p - 2 ) S n>0 )] 

oo m 
s=0 1=0 

+ 2ir 2 fl r (k)5 nfi (5 Pjl - <J Pi _i) . 
The following constants were introduced: 

-Is 

dr^-j Tji [i"m\J ji) '° ^v{k Tji) j 



H-rn,u{k) 
^a,u{k) 

fx r (k) 

r"m,a 
An 

-pn 
1 si 

Otfi, 



»(2-+l)(jfe) 



dr^i rjj ^t a (rji) J u{k fji) 



^ r ji r ji ^riXjijJlik Tji) 
dr^j rji flm,a(j'ji) \k fji) 



2s+l 



2 2s (s -l)\(s + 1 + 1)1 

:-i) a i 

2 2s (s -Z)!(s + / + 1)! 

2tt 



«n-£-l 



(2/-l)(2/ + l) (2/ + 3)(2/ + l) 

&n+l Ctn+l+1 



d£ sin I - ) e" ir * = 

1 2 / 1 - 4n 2 



(2/-l)(2/ + l) (2/ + 3)(2/ + l) 
4 



[n e Z). 



(A.6) 

(A.7a) 
(A.76) 
(A.7c) 
(A.7d) 
(A.7e) 

(A.7j) 
(A.7^) 



The Bessel functions of the first kind are denoted by J u (x). Please note, that (1A.6j) is 
only valid for k = (k, 0). 



Appendix B. Transport Coefficients of (14. lift 

The transport coefficients occuring in the hydrodynamic theory (14. lip are listed below. 



Ai = CiP - £5 

2" \ ~9~ P ~ P ~ S ° 



(B.lo) 
(B.16) 
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